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Alternative statistical-mechanical descriptions of decaying two-dimensional turbulence 

in terms of "patches" and "points" 
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Numerical and analytical studies of decaying, two-dimensional (2D) Navier-Stokes (NS) turbu- 
lence at high Reynolds numbers are reported. The effort is to determine computable distinctions 
between two different formulations of maximum entropy predictions for the decayed, late-time state. 
Though these predictions may be thought to apply only to the ideal Euler equations, there have 
been surprising and imperfectly-understood correspondences between the long-time computations 
of decaying states of NS flows and the results of the maximum-entropy analyses. Both formulations 
define an entropy through a somewhat ad hoc discretization of vorticity to the "particles" of which 
statistical mechanical methods are employed to define an entropy, before passing to a mean-field 
limit. In one case, the particles are delta-function parallel "line" vortices ("points" in two dimen- 
sions), and in the other, they are finite-area, mutually-exclusive convected "patches" of vorticity 
which in the limit of zero area become "points." The former are taken to obey Boltzmann statis- 
tics and the latter, Lynden-Bell statistics. Clearly, there is no unique way to reach a continuous 
different iable vorticity distribution as a mean-field limit by either method. The simplest method of 
taking equal-strength points and equal-strength, equal-area patches is chosen here, no reason being 
apparent for attempting anything more complicated. In both cases, a non-linear partial differen- 
tial equation results for the stream function of the "most probable" or maximum-entropy state, 
compatible with conserved total energy and positive and negative vorticity fluxes. These amount 
to generalizations of the "sinh-Poisson" equation which has become familiar from the "point" for- 
mulation. They have many solutions, which must be obtained numerically, and only one of which 
maximizes the entropy on the basis of which it was derived. These predictions can differ for the 
point and patch discretizations. The intent here is to use time-dependent, spectral-method direct 
numerical simulation of the Navier-Stokes equations to see if initial conditions which should relax 
to different late-time states under the two formulations actually do so. 

PACS numbers: M 47.10. +g. M 47.11. +j. M 47.27.Jv. M 47.32. Cc. 



I. INTRODUCTION 

It has been known for several years that two- 
dimensional Navier-Stokes (2D NS) turbulence with peri- 
odic boundary conditions and at high Reynolds numbers 
(in excess of a very few thousand) will relax to long- 
lived, quasi-steady states whose topology is preserved 
and whose energy decays more or less inversely propor- 
tionally to the Reynolds number, computed with respect 
to a box dimension and an rms initial turbulent velocity. 
When this energy decay time is large compared to the 
initial eddy turnover time, the quasi-steady state can be 
reached at a time when the total energy decay is fraction- 
ally small, though the enstrophy decay can be fractionally 
large. It came as a surprise when, a decade or more ago, 
a series of such computations ||, ||, [| revealed that the 
late-time quasi-steady state for an initially turbulent run 
at a Reynolds number above 14,000 showed a pointwise 
hyperbolic sinusoidal dependence between stream func- 
tion and vorticity. 

The prediction of such a dependence, in the con- 
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text of a mean-field treatment of ideal line vortices 
(or guiding-center plasma rods) had been given thirty 
years ago [| | and has since been extended and 
refined in a series of investi gati ons by several groups 
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every case referring to ideal, non-viscous systems. The 
system is Hamiltonian with a finite phase space, and it is 
natural to apply Boltzmann statistics to its dynamics, as 
originally suggested by Onsager (see also Lin p2[). 
The surprise came in the extent to which the ideal Euler 
mean-field predictions fit the Navier-Stokes results. At 
least one attempt was made to define an entropy for the 
case of finite viscosity |2j| but generated new puzzles of 
its own. 

In the late 1980s and early 1990s, an alternative for- 
mulation was given by Robert, Sommeria and Chavanis 
||, [27), by Miller and colleagues Jpj and later 
explored by Brands, Maassen and Clercx [|29fl. The prin- 
cipal difference was that the vorticity field was discretized 
not in terms of delta-functions, but rather in terms of 
finite-area, mutually-exclusive "patches" of vorticity, to 
which Lynden-Bell statistics |jo| could be applied. The 
choice of parameters in the patch formulation is wider 
than that of equal-strength point vortices, in that one 
must decide in advance the size of the patches, and the 
number of "levels," or strengths, that the patches carry. 
There seems to be no deductive mechanism for making 
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such choices, and so we stick with the simplest, taking 
equal-area patches and never more than three levels, in- 
cluding zero. 

For readers not inclined to pursue the original refer- 
ences for these formulations, which in some cases require 
considerable mathematical sophistication, we offer in Ap- 
pendix A a summary treatment of the derivation of the 
"patch" dependence of vorticity on stream function in 
the "most probable," or maximum-entropy, state. It is 
attempted to do this using the least forbidding and most 
transparent mathematics available to the purpose. We 
show how in the limit of zero- area patches, the point de- 
pendence is recovered. Thus the reader is being referred 
to Appendix A for a derivation of the nonlinear partial 
differential equations to be introduced in Section |n| as 
candidates for the predictor of the most-probable state 
of late-time decaying 2D NS turbulence. 

Section || will summarize both the numerical method 
(now well-known in its essentials), due originally to 
Orszag and Patterson |3ll [32] , to which various sets of 
initial conditions will be subjected. We will also describe 
the motivation for choosing these initial conditions, in 
terms of what we expect from them as possible most- 
probable states. We will concentrate on initial conditions 
for which it appears that the predicted outcomes using 
point and patch entropies will be as different as possible. 
We will refer the reader to Appendix B for a description 
of how the sinh-Poisson relation and its patch general- 
izations are dealt with numerically (a non-trivial task). 



Section 111 contains the main body of the results that 
we have found. Among other unexpected features, we 
have found a series of one-dimensional (ID) solutions to 
the most-probable state equations, whereas the literature 
up to now has dealt only with 2D solutions ("dipoles," 
"quadrupoles," "octupoles," etc.). There are cases in 
which the entropies of the ID solutions (which we call 
"bars" ) are extremely close to those of the dipole over a 
considerable range of energy. Which is greater depends 
on seemingly arbitrary choices, in the patch formulation, 
such as the areas of the patches used. Classes of initial 
conditions are found (not essentially turbulent ones) for 
which relaxation to the "bar" states is observed. It will 
be easier, in the context of the details of the solutions, 
to illustrate the variety of behavior we have been able to 
catalogue. 



Section IV will be a summary of what we think we have 
learned, and of what remains to be learned. One of the 
more radical, if incidental, conclusions which we believe 
may be extracted from these runs (to be discussed in 
more detail later) is that in the context of the initial value 
problem, all 2D NS turbulence may result only from its 
having been there initially. Otherwise, it can come into 
existence only through random external "stirring." 

We should caution the reader that the word "turbu- 
lence" in this manuscript will be used in a somewhat 
looser sense than is customary. It will not always mean 
a state of the fluid in which the kinetic energy is widely 
shared by many Fourier modes. We have concentrated 



on initial conditions which seemed to us most relevant 
to producing effects peculiar to the "patch" description: 
namely, initial vorticity distributions in which a few flat 
levels of nearly constant vorticity could be readily iden- 
tified. We expected, because of instabilities in the shear 
flows these represent, that turbulence of a typical broad- 
band modal structure would soon be generated, as it is 
in three dimensions. We found that instabilities could 
indeed be activated, but often, broad-band turbulence of 
the conventional kind was not the result. Instead, energy 
spectra dominated by only a few low-lying (in Fourier 
space) modes were nonlinearly dynamically converted to 
Fourier spectra dominated by only a few (but different) 
low- lying modes. We have in fact found it nearly im- 
possible to generate genuinely broad-band 2D NS tur- 
bulence through instabilities, rather than through initial 
conditions or externally-imposed stirring. Perhaps sur- 
prisingly, the dynamical evolution involved sometimes, 
but not nearly always, led to late-time quasi-steady states 
that seemed to correspond to most-probable "patch" pre- 
dictions: most interestingly, the one-dimensional "bar" 
state, as will be seen. Broad-band, initially-excited tur- 
bulence continued to lead to the classical dipolar late- 
time state, as it has in previous turbulence simulations. 



II. THE NUMERICAL AGENDA 

We address ourselves to the long-time dynamics in two 
dimensions (2D) of the Navier-Stokes equation in the 
usual vorticity representation, 



dui 
~dt 



zA7 ui, 



(1) 



where the fluid velocity v is to be written in terms of 
the stream function ip as v = V x (e-0), and e is a unit 
vector in the z direction, v has only x and y components, 
which depend only upon x, y, and the time t. In the 
natural dimensionless units of the problem, the kinematic 
viscosity v may be interpreted as the reciprocal of the 
Reynolds number, which we will specify in more detail 
presently. The curl of v is the vorticity oj = (0,0, ui), 
which has only one component, also in the z-direction. 
The stream function ip and the non- vanishing component 
of vorticity lo are related by Poisson's equation, 



vV = 



(2) 



We note that dropping the final viscous term in Eq. (jl ) 
leaves us with the 2D equations of an ideal Euler fluic . 
We note moreover that a time-independent solution of 
these Euler equations results any time uj is a differentiable 
function of ip. 

We will work in a periodic box with sides 27r in length, 
and will use the unit length 1 to define the Reynolds 
number; thus our basic unit of length is roughly 1/6 of 
a box dimension. The unit of velocity will typically be 
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a root-mean-square value of the initial velocity, and we 
will attempt to make this equal to 1 whenever possible. 
Thus v in Eq. ([[]) may reasonably be identified with the 
reciprocal of an initial Reynolds number, Re. We are 
interested in values of Re of at least several thousand, 
so that the final term in Eq. (0) is formally very small, 
except in regions of steep vorticity gradients. The "eddy 
turnover time," which we shall use as a unit of t, is thus 
about 1/6 of a box dimension divided by an initial rms 
velocity. It has been known for many years that under 
such circumstances, the kinetic energy E, defined by 



E = - 
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decays slowly proportionally to v. However, the decays 
of higher-order Euler-equation invariants constructed as 
integral moments of u>, such as the enstrophy, 



n = i 



i 



2 (2tt) 2 
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appear to continue to decay at an 0(1) rate; even a weak 
dependence of their decay rate upon v has not been con- 
vincingly demonstrated. Note that both ideal invariants 
are referred to unit volume, in the dimensionless units. 

When high Reynolds numbers exist, then, the turbu- 
lent decay of energy is seen to be very slow, relative to 
any other identifiable ideal Euler invariant. When the 
energy decay time is large enough to be well-separated 
from the eddy turnover time, it is possible to observe 
in numerical computations [jl], |2|, ^ that a quasi-steady 
state is reached in a time over which the fractional decay 
of energy is small (a few hundred eddy turnover times). 
One might expect to find a "selectively decayed" state, 
in which the enstrophy to energy ratio is minimal p3[ , 
and indeed, if one waits long enough, it can be analyti- 
cally proved that such a state must be approached [p4[ . 
However, it was found || that long before that time, 
the quasi-steady, slowly-decaying state that is reached 
has a rather sharp one-to-one pointwise correspondence 
between u> and if>. The elucidation and testing of this 
correspondence is the principal purpose of this paper. 

Clearly, even though a slow viscous decay may be su- 
perimposed on such a state with u> approximately a func- 
tion of ip, the state is also closely approximated by a 
time-independent solution of the Euler equations. There 
is no a priori reason why this should be so. It is a fact 
we attempt to incorporate here into a coordinated set of 
direct numerical simulations of Eq. ([!]) and a combined 
analytical and numerical argument, based upon statisti- 
cal mechanics, in pursuit of what the connection between 
uj and tp should be. 

The reader is referred to Appendix A for a summary of 
the statistical mechanical arguments. They depend upon 
a discretization of vorticity in terms of delta-function 
line vortices ("points") or in terms of mutually-exclusive, 
finite-area "patches." Boltzmann statistics are applied 



to the former, to define an entropy, or logarithm of the 
probability of a state, and Lynden-Bell statistics are ap- 
plied to the latter. The limit of zero-area, finite-vorticity 
patches are points. Thus the patch formulation can be 
viewed as containing the point version of the theory, 
and either is or is not a useful generalization of it. In 
both cases, what results is a "most probable" dependence 
of vorticity upon stream function. A mean-field limit 
(infinitely-many points or patches, of arbitrarily small 
strength) is then taken, to yield continuous and differen- 
tiable functions co(ip). In Appendix A, it is shown that 
for points, the resulting function is: 
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For patches, the resulting function is: 
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1=0 



(5) 
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The symbols a, (3 are Lagrange multipliers, which en- 
ter via a maximization of the appropriate entropy subject 
to the constraints of given energy and positive and nega- 
tive fluxes of vorticity. They are determined in principle 
by demanding that the energies and vorticity fluxes cal- 
culated on the basis of Eqs. (||) and (||) or (g) match 
specified values. In Eq. (f3|), the Kj are the "levels" of 
vorticity corresponding to the different sized patches, and 
must be chosen somewhat arbitrarily. A/M is a fixed size 
of a "patch," defined in detail in Appendix A, and may 
be chosen arbitrarily, within wide limits. 

Both Eqs. (||) and (^J) have infinitely many solutions. 
The physical ones are interpreted to be those which max- 
imize the appropriate entropy from which they were de- 
rived. The others represent local maxima, which, as we 
shall see, may in some cases represent attainable states 
in the computations. 

The questions before us are: 

• How does one determine the entropy-maximizing 
time- independent solutions to Eqs. (||) and (g)? 

• How close may we come to those solutions in a 
dynamical computation that solves Eqs. (Q) and 
(0) as a time-dependent initial-value problem with 
large but finite Reynolds number? 

• Are there noticeable differences between the pre- 
dictions of Eqs. (|B|) and (||) and are there reasons 
for preferring one to the other as predictors of late- 
time turbulent decays? 

It is to be stressed that all three questions are answer- 
able not in the abstract, but only as a result of somewhat 
demanding computations. There is no a priori reason 
why Navier-Stokes turbulent decays should be predicted 
at all by anything having to do with the Euler equations. 
The latter will never be soluble, for continuous initial 
conditions, over time intervals long enough to make ideal 
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solutions of much interest, because it is in the nature of 
Euler codes, having no minimum physically-determined 
length scale, to overrun their own resolution in a very 
few eddy turnover times. If the predicted states had not 
shown some empirical relevance to Navier- Stokes solu- 
tions, there would be no justification for this activity. 

The dynamical code used here is of the now familiar 
Orszag-Patterson pseudospectral variety 
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The 

code is a parallelized MPI Fortran 90 version of an earlier 
Fortran 77 code provided by W.H. Matthaeus (private 
communication). It is fully de- aliased, using the shifted- 
grid method. It has been run, in the runs reported here, 
on the SGI Origin 3800 at SARA Supercomputing Cen- 
tre in Amsterdam. Some of the simulations (resulting in 
the "bar" final state in section III B 1 b| ) have also been 
recomputed with a different pseudospectral Fourier code 
(provided by Dr. A.H. Nielsen, Riso National Labora- 
tories, Denmark.) These computations yielded the same 
conclusions. 

All runs resolve the Kolmogorov dissipation wave num- 
ber based on enstrophy dissipation. Previously, ini- 
tial conditions for turbulent decay runs have tended to 
use randomly loaded Fourier coefficients in the spec- 
trum or vorticity fields up to some upper wave number. 
These have typically been chosen to match some cascade 
wavenumber fc-spectrum, such as the "—3" direct enstro- 
phy cascade |3^] spectrum predicted by Kraichnan 
p5[ . The phases have been chosen from a random num- 
ber generator. This has been done to achieve the most 
disordered (in some sense) initial field compatible with a 
particular power spectrum. Several of the runs to be re- 
ported later are in fact of this variety. Such an initializa- 
tion can only yield a vorticity distribution that is analytic 
in x and y, since only sinusoidal functions are involved, 
even though the spatial dependence may be wildly fluctu- 
ating. Such a vorticity distribution can take on any par- 
ticular value of ui only over a set of measure zero. For this 
reason, one might imagine it would tend to de-emphasize 
features that corresponded to the "levels" in the patch 
formulation, over which the vorticity is supposed to take 
on a constant value inside a compact area. Therefore, we 
have also stressed initial conditions that have large, flat 
areas of vorticity, separated by thin regions, with as steep 
spatial gradients between them as the code will resolve. 
Our original intention was that these would correspond 
to unstable laminar shear flows whose instabilities would 
subsequently generate turbulence. Somewhat to our sur- 
prise, we have found it easy to produce shear-flow insta- 
bilities but not ones that led to turbulence, in the sense 
of broad fc-spectra. We have come to suspect that this is 
a feature inherent to two-dimensional flows, and perusal 
of the literature has uncovered only three-dimensional 
turbulence as a consequence of unstable shear flows, but 
not two-dimensional; but this must be left as a conjecture 
rather than a demonstrated fact. 

In any case, in order to induce a large number of 
Fourier modes to participate in the subsequent dynam- 
ics, we have found it necessary, in the runs evolving from 



vorticity distributions with flat areas, to add significant 
amounts of random noise initially, in order to get a sub- 
sequent fc-spectrum that could readily be called "turbu- 
lent." Even then, some of the evolution we find might 
have its status as turbulence disputed. 

So referring to Appendix A for the derivations of Eqs. 
(||) and (|J), and to Appendix B for the method of their 
numerical solution, we proceed in Sec. 3 to a description 
of initial conditions used in the dynamical runs, the states 
to which they evolved, and their comparisons with the 
maximum-entropy or "most probable" states that Eqs. 
(H) and M) imply. 



III. NUMERICAL RESULTS 

A. Extracting Maximum-Entropy ("Most 
Probable") States 

We first summarize the results of the solutions of the 
most-probable-state equations, derived in Appendix A; 
the method of their solution is described in Appendix 
B. The goal is not only to find numerical solutions but 
to determine which among the solutions has the highest 
entropy for given energy and vorticity fluxes. Generally 
speaking, the entropies of states with more maxima and 
minima have been seen to be less than those with fewer, 
and we can be guided by this. But there are instances 
where solutions of two different topologies can be found 
whose entropies lie very close to each other, as found by 
Pointin and Lundgren j|] , and we must concentrate on 
those. 

We consider the point relation, Eq. (|B|) in the absence 
of any conditions that would suggest asymmetry between 
positives and negatives, so that the two Lagrange mul- 
tipliers a + = ot- = a. The result is the sinh-Poisson 
equatio n (a typical lu — ijj dependence of which is shown 

in Fig. nm 



= 2e Q sinh(/?V), P < 0. 



(7) 



We also consider the patch relation, Eq. (^), special- 
ized to the three-level case of vorticity levels —1, 0, and 
+1. Again ass umin g symmetry between positives and 
negatives (Fig. 1(b) ), we have 
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where (3 < and A /M is the (arbitrary) size of a patch. 
Eq. (m) can be rewritten as 



V 2 * = -A 2 sinh 4- 
where we have defined 
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FIG. 1: Typical u) — if) relation from (a) sinh-Poisson, (b) 3-level Poisson equation and (c) tanh-Poisson equation (special case 
when e~ a -> in Eq. (|)). 
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FIG. 2: Contours of constant stream function tp f° r the solutions of Eq. (j5|). Negative values are shown as broken lines, 
throughout. 



and 



A = 2 I 



(11) 



The form of Eq. (Q) makes it most susceptible to nu- 
merical solution. Once it is solved, for a given value of 
A, we must keep in mind that in order to evaluate the 
entropy associated with the solution, we must revert to 
the parameters of Eq. ([?]), while guaranteeing that a and 
/3 satisfy Eq. (|Tl|). Since for each fixed value of A 2 , we 
have infinite sets of possibilities for a and /3, a recipe is 
needed for choosing them. Since our goal is to plot the 
entropy vs. energy for a fixed value of the vorticity flux, 



-Mdxdy 



for the domain [0, L] x [0, L] (here we let L = 2n and 
VL$ — 1, for convenience), we must also be assured that 
the solution we get for ip satisfies this condition. 

Fortunately, these two conditions can be satisfied si- 
multaneously. Combining Eqs. (|l0|), (11) and = 1 
leads to the conclusion that 1/31 = 
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Because Eq. (||) is the equation solved, the value of 
il^ is readily obtained. Thus we obtain /3 and from Eq. 
(|lT|), a. We have now all the parameters that correspond 
to fixed vorticity flux of either sign. From Eqs. ( |l3| ) and 
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by letting: 
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we can draw plots of the entropy vs. energy (E^) 
for fixed unit flux of positive and negative vorticity. 

It is useful to introduce terms for the solutions of dif- 
fering topology. For example, for a given energy and vor- 
ticity flux, the cont ours of stream function may look like 
Figs. |2| Fig. 2(a) will be called the "dipole" solution. 
Fig. 2(b)| will be called the "quadrupole" solution, and so 
on through higher "multipoles" that correspond to suc- 
cessively higher (even) numbers of maxima and minima. 
We have also found one-dimensional solutions as illus- 



trated in Fig. 2(c). These one-dimensional solutions will 
be called "bar" solutions, and there are also an infinite 
sequence of them, with basic periodicities 2ir, tt, 7r/2, ... 
These states are obtained, as explained in Appendix B, 
by iterating a trial solution with similar topology until 
convergence is obtained. 




FIG. 3: Entropy vs. energy for the solutions of Eq. (JsJ) for 
unit positive and negative vorticity flux, computed for the 
"point" discretization. 

Fig. H shows a plot of the entropy vs. energy, for 
fixed unit positive and negative vorticity flux, for the 
quadrupole, bar, and dipole solutions. Observe the very 
small difference that makes the entropy of the dipole 
greater than that of the bar. It will turn out that it 
is possible to find either solution as a consequence of the 
development of certain initial conditions, as we shall see 
later. 

Turning to the patch prediction, Eq. (|^) can be sim- 
plified as 
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Solving Eq. ( |l5| ) is not as easy as solving Eq. (||) 
because there is an extra parameter, g, which implictly 
reflects the arbitrary choice of the patch size. Another 
strategy is required: 

1. We choose the size of the patch, so that 
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2. Then a trial value of g is chosen before solving Eq. 
(p"5|), again iterating about a trial solution of desired 
topology. Using Eq. (|l^), we have the value of (3, 
and a from Eq. (|l^), and Q^, are given by 
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(21) 



However, the parameters obtained in this way are 
not usable since the condition of unit positive and 
negative vorticity flux is not in general satisfied. 
We must return to the second step and change g 
until the desired accuracy is obtained from Eq. ( |2l| ) 
1. Once this has been done, we are 
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then in a position to evaluate the entropy from the 
algebraic expression (see Appendix A): 
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We may plot the entropy of Eq. (J22J) in Fig. 4(a), 
which is obtained by the choice of M/A = 3.7814, a 
relatively large "patch" size. It will be seen that for this 
large patch size, the entropy of the bar solution is greater 
than that of the dipole, a different conclusion than that 
of the point calculation. However, if we reduce the size of 
the patch, we find that the result approaches that of the 
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point formulation. In Fig 
M/A = 25, and in Fig. ||c 



4(b), we see the result when 
the result when M/A = 100, 
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where the maximum entropy state is again the dipole one. 
Considering the very small entropy differences between 
the dipole and the bar solutions shown in Figs. || and 
U, it might not be thought surprising if the fluid had 
difficulty making up its mind which state to relax into. 

Before turning to the results of the dynamical compu- 
tations, we offer a few observations on the relations of 
the patch versions of the theory to the point version, and 
to each other. It is clear from the derivation that in the 
limit that the patch sizes become smaller and smaller, at 
fixed and finite separation, the point version of the theory 
is recovered. Depending upon the number of levels cho- 
sen for the patch formulation, there are many versions of 
the most-probable-state patch equation; the three-level 
version of Eq. (|§|) is not the most general by any means. 
Each will have different solutions. There is, however, no 
apparent unique or practical prescription for how many 
levels, or what size patches, should be used to represent a 
particular initial analytic vorticity distribution with high 
accuracy. It is not clear that it can be done without re- 
quiring the patch size to shrink to zero, at which point 
it becomes indistinguishable from a point representation. 
In fairness, we should also say that there is another aspect 
to the point formulation that is also ambiguous: namely, 
there is no reason to choose the point vortices of the 
mean-field theory to be of equal strength. If some are 
of different strength, Eq. (Q) will also change. Keeping 
these in mind, we turn now to the effort to see some of 
the solutions as consequences of direct numerical solution 
of the Navier-Stokes equation. 
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FIG. 4: Entropy vs. energy at unit vorticity fluxes for the 
solutions of Eq. (^): (a) with a relatively large patch size 
(M/A = 3.7814); (b) with a somewhat smaller patch size 
(M/A = 25); (c) with a still smaller patch size (M/A = 100). 
Note that in (c) the dipole solution has become slightly more 
probable than the bar. 



B. Dynamical Solutions and Comparisons 

We now present the results of dynamical, pseudo- 
spectral-method solutions of the 2D NS equation, using 
a resolution of 512 2 . The time step is 0.0005. The initial 
energy, using the normalization of Eq. (^|), is 0.5. There 
is no hyperviscosity or small-scale smoothing of any kind. 



1. "Dipole" and "bar" final states. 

a. Random, broad-band initial conditions. The first 
run we report here is essentially a reproduction of the run 
of Matthaeus et al. jl], ||, || , but with a lower Reynolds 
number (1/v = 5000), corresponding to an initial Taylor- 
scale Reynolds number 




558. 



R\ has increased to 3143 by the end of the run. 

This run served as a benchmark for our parallelized 
code in the beginning, and the McWilliams initial condi- 
tions fl38| can also be used as noise, with variable ampli- 
tude, to be added to later simulations to break unwanted 
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FIG. 5: Equally spaced contours of constant vorticity at three 
different times (left column) and corresponding modal ener- 
gies at the lower values of k (right column) , during the evolu- 
tion of the McWilliams/Matthaeus initial conditions. These 
have no flat patches of vorticity initially, even approximately. 



symmetries and accelerate the dynamical development. 
The run also provides an opportunity to introduce the 
types of data displays that will be used throughout the 
rest of this Section. The Fourier modal energy spectrum 
E(k) = (1/2) |v(k)| 2 is initialized according to 



E(k) 



C 



1 + (|) 4 



for 1 < k < 120 (here, and for this purpose only, the wave 
number k is binned in integer values by a standard For- 
tran routine, and the partition of energies among modes 
with the same integer k is decided by a random number 
generator), and zero otherwise. The phases of the Fourier 
coefficients are chosen from a Gaussian random number 




-2 -1.5 -1 -0.5 0.5 1 1.5 2 

FIG. 6: The lo — ip sc atter plot for the run shown in Fig. || 
(which is close to Fig. 1(a)). 



generator. C is a constant to be adjusted to make the 
total initial energy equal to 0.5. 

The results were quite similar in all respects to those of 
Matthaeus et al. jlj [| [|, so only a few figures (Figs. || 
and ^ are shown here. The contour plots of vorticity and 
stream function relax to the familiar dipole final states. 
The left column of Figs. || shows the vorticity contours at 
three separate times. The right column of Figs. [5] shows 
the modal energy spectra, for the lower part of fc-space 
(the maximum k is 241) at those same times, exhibited 
as three-dimensional perspective plots. They are initially 
broad-band, but evolve as expected to be concentrated 
in the lowest values of k. In Fig. [| a scatter plot of u> vs. 
if) shows a good agreement with the hyperbolic sinusoidal 
dependence of vorticity on stream function, as predicted 
in the point formulation (Fig. 

b. "Bar" A different, and interesting, evolution re- 
sults = 5000) if we initialize using the quadrupole 
solution from Eq. ([l5"|). In this run, R\ increases from 
2391 initially to 4920 at the end. This is not predicted 
to be the maximum entropy state from Eq. (|2^) , or any 
other equation, patch or point; and so should evolve when 
initialized with some random noise. The quadrupole so- 
lution was found to have a symmetry which persists in 
time, and the noise, it is hoped, will permit that sym- 
metry to be broken. It was found that the symmetry 
persisted until the noise was raised to quite substantial 
levels. Figs. fj] are perspective plots of the initial vortic- 
ity as a function of x and y, with and without the noise 
added. The noise level in this case is about twice the 
quadrupole vorticity, and its randomness should break 
any symmetry that might be present. 

Fig. H shows the evolution of the vorticity and stream 
function contours that result from the initial conditions 
shown in Fig. |7(b)| . The intermediate panel (t = 15) 
shows no obvious connection to the symmetry of Fig. 
7(a) or to that of the final panels (t=1250) of Figs. ^ 
(the first is odd about the center lines of the basic square, 
but clearly, the t=15 state has no such symmetry). The 
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(a) 



FIG. 7: The initial vorticity field as a function of x and y for 
a run intended to exhibit patch characteristics, (a) without 
random noise, (b) has had substantial random noise added to 
the vorticity field. 



relaxation to the one-dimensional "bar" state at t=1250 
has been quite robust in several such runs. Fig. ^| is 
a ui — ip scatter plot for the bar state. There has been 
considerable decay of vorticity, as evidenced by the dis- 
appearance in Fig. ^ of the high and low values visible 
in Fig. |7(b)| . By t = 1250, the energy (E) has decreased 
to 58.4% of the initial value. 

When we consider the modal energy spectra in Fig. [l^, 
for this evolution, we see that both the initial and final 
states are dominated by four and two Fourier modes, re- 
spectively. In fact, the evolution after t=40 has this char- 
acter. The energy spectrum at t=15 is somewhat more 
broad-band than either, but whether it should properly 
be called "turbulent" may be debated; it never achieves 
the fully broad-band character of the evolution shown in 
Fig. ||, for example. The final state achieved is consistent 
with the maximum-entropy prediction of the " patch " for- 
mulation for suitably large-sized patches (Fig. 4(a) ) , but 
not w ith the same formulation for smaller-sized ones (Fig. 
f4(c)| ). The right column of Figs. |l^ shows a histogram 
of the vorticity at three different times in this evolution; 
it is in effect a frequency distribution for the appear- 
ance of various values of the vorticity over the plane, one 



per computational cell in the 512 array. It will be seen 
that this distribution changes enormously over the course 
of the run, which it cannot do, of course, in any ideal 
Euler picture. Figs, [ll] are binned plots of the angle- 
averaged modal energy spectrum (E(k)) and enstrophy 
(fi(fc)) spectrum at the three different times. Fig. |l2| 
shows the time history of several global quantities for 
this run: energy, enstrophy, mean wave number (square 
root of the ratio of enstrophy to energy) , and square root 
of the ratio of palinstrophy (P = E k (l/2)fc 2 |w(k)| 2 ) to 
enstrophy, showing spikes that are characteristic of vor- 
tex merger events at high but finite Reynolds numbers. 

Figs. |7(b)| through [lj, show, then, a reproducible 
evolution of an initial condition consisting of a patch 
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FIG. 8: Contours of constant vorticity (left column) and con- 
stant stream function (right column) at three different times 
for the run originating fro m th e intial vorticity distribution 
(with noise) shown in Fig. 7(b). 
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FIG. 9: The u — ip scatter plot of for the late-time state 
achieved in the run shown in Fig. @ (which is close to Fig. 
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FIG. 10: Low-k modal energy spectra (left column) and vor- 
ticity histogram (ri ght c olumn) at three different times for the 
run shown in Figs. 7(b) and ^. The vorticity histogram shows 
the number of times a particular value of uj is recorded as one 
cycles through the computational grid. Notice that there is 
substantial variation from one time to another. (An Euler 
equation solution should preserve this histogram in time.) 





FIG. 12: Time evolution, for the run shown in Figs. 7(b) 
through |n], of the four global quantities energy, enstrophy, 
\JQ/E or mean wave number, and ratio of the square root of 
palinstrophy to enstrophy. Note that the mean wave number 
(Taylor microscale) is never large, and so the strict defini- 
tion of turbulence is not fulfilled at any time during the run. 
Nevertheless, the topology changes dramatically. This is ef- 
fectively a selectively decayed state J32j ], and may not invite 
a probabilitistic explanation. 



quadrupole plus large amounts of random noise into a 
one-dimensional "bar" most-probable state. This state 
is more probable than the dipole, according to a "patch" 
analysis if the patch size is big enough, though it is not 
for a smaller but finite patch size. 

A similiar initial condition is used by Segre and Kida 
p6| . However, their computation made use of hypervis- 
cosity and appears not to have run long enough for the 
proper late-time state to evolve. 
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FIG. 11: Angle- averaged modal energy and enstrophy spectra 
at three different times for the run shown in Figs. H, Pi and 
[Ifj| . Note that these are not the usual omnidirectional spectra; 
since they do not contain the factor 2irk, their maxima occur 
at the minimum k. 



2. Local maximum states 

a. Evolution from a 64~pole initial condition. The 
contrasting evolutions of the two initial states in the 
previous sub-section naturally arouse curiosity about 
whether there is a limit in which one behavior goes over 
into the other. In this sub-section, we attempt an answer 
to this by considering initial conditions which originate in 
high-order multipole solutions of the patch formulation, 
not maximum entropy states in either formulation, but 
intuitively closer to the random initial conditions that led 
to the dipole solution before, in the first run reported. 

There are four numerical solutions in this group, with 
1/V = 10, 000 for all four runs. They all originate in the 
64-pole solution of the patch formulation, using the 3- 
level equation, Eq. (|l5|). The same conclusions that we 
reach can also be reached using 16-pole initial states, but 
we will not display those results here. 

The motivation was to see if the high-order multipole 
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FIG. 13: Time evolution of a 64-pole initial condition, triggered only by round-off error, into a bar state (the scatter plot of 
final state is close to Fig. 1(c)). No noise beyond round-off error has been added to the initial condition. 
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FIG. 14: Time evolution of the vorticity contours, starting from the same initial condition as in Fig. |13| but with a healthy 
addition of random noise (In this run, R\ increases from 2036 i nitiall y to 11400 at the end.). The last panel is the late-time 
w — tj) scatter plot for the resulting dipole (which is close to Fig. 1(a) ). 
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FIG. 15: Evolution of the one-dimensional "8-bar" initial condition, with random noise added initially (In this run, R\ increases 
from 3228 initially to 11500 at the end.). The evolution is toward a dipole, now, as seen in the bottom panel (which is close to 

Fig. ra. 



initializations would lead to "bar" final states, the way 
the quadrupole does. Intuitively, they would seem to be 
closer to our picture of what true turbulence might look 
like. The bar states are no longer found, but there are 
"local maximum" entropy states that can be achieved in 
the limit of low noise in the initial conditions. 



First consider what happens to the 64-pole initial con- 
ditions without any noise, as shown in Fig. [l3| A 
straightforward laminar evolution occurs, with the end 
product (t = 250) being a one-dimensional bar state with 
a total of eight maxima and minima. This state is essen- 
tially dominated by one Fourier mode, as indicated by 
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FIG. 16: Evolution from initial conditi ons c orresponding to a sinh-Poisson quadrupole, plus random noise. Note the differences 
in evolution from those shown in Figs. 7(b) - |l^, whose initial conditions are superficially similar. The low-k part of the Fourier 
space contains most of the energy throughout, so the evolution again is not classically "turbulent." 
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FIG. 17: Vorticity histogram at three different times for the evolution shown in Figs. |l6|. 
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FIG. 20: Evolution of a 4-bar initial condition, plus random noise, that ceases to evolve (perhaps because of loss of nonlinearity 
due to decay) before either a dipole or a bar final state is achieved. 



the essentially linear pointwise dependence of vorticity 
on stream function. In Figs. [l4|, we show the evolution 
of the same initial conditions with only a low level of 
noise: 1/2000 of the amplitude of the 64-pole solution, 
not visible on the t — contour plot. In this case, fully- 
developed turbulence does seem to result, with a dipole 
final state as the end result. 

The third run in the group concerns the result of 
taking the final bar state (we raise the initial energy: 
E(t = 0) = 0.5), without putting any noise into it, from 
Figs. 
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and allowing it to run. This appears to be a 
time-independent state, stable in the presence of round- 
off error for the duration of the run, but not a maximum- 
entropy state and not stable in the presence of noise of 
greater amplitude. This is clear from Figs. |l5|, which 
show the development of the unstable evolution induced 
by putting on the same level of random noise as in Figs. 
|l4| , with a fully-developed turbulence and a dipole final 
state as the result. 

There are thus some subtleties revealed by these four 
runs. Bar final states, predicted by the patch theory with 
sufficiently large patch size, do result from the evolution 
of a patch quadrupole plus sufficient noise. However, 
higher order multipoles from the patch formulation seem 
to be unstable at low levels of random noise, and evolve 
into dipoles, consistently. Still less easy to fit into the 
picture is the laminar evolution of the "local maximum" 
64-pole state without noise into the bar state shown in 
Figs. |l^. Neither state is in this case an absolute max- 
imum entropy one, though both are in some sense local 
maxima. 

b. Evolution starting from a sinh-Poisson quadrupole. 
We have started several runs from initial conditions that 
originate with a sinh-Poisson quadrupole state, with 
\jv = 10,000. The first is the quadrupole state without 
any added random noise, and the others are the same 
state plus smaller or larger amounts of random noise. 

The zero-random noise initial conditions seems to be 
stable and to remain in the same shape in the presence 
of only round-off error for the duration of the run, as 



reported previously J37j|. Even after putting in enough 
random noise to manifestly break the symmetry (Figs. 
|l?j| ) , the quadrupolar shape is maintained for quite a long 
time before breaking down into a dipole by t = 300. Figs. 
|l7| are vorticity histograms computed at three different 
times, showing more concentration of uj near zero values, 
characteristic of sinh-Poisson states. The oj — ip scat- 
ter plots shown in Figs. [l8| also illustrate this feature. 
Evolution of the global quantities for this run is shown 
in Figs. |l9], showing abrupt changes in enstrophy and 
palinstrophy when vortex merger occurs, but otherwise 
not displaying characteristics of turbulent behavior. Note 
that energy is well conserved for this case. In this run, 
the quadrupole persists for a long time because the en- 
ergy is concentrated in the four well-separated vortices, 
before breaking down into a dipole. This is noticeably 
different behavior than when we started with the patch 
quadrupole. 



3. Oddities 

a. A run (\jv — 20,000) that leads to an "unclas- 
sifiable" final state is shown in Figs. [2(]. It begins with 
a four-fold bar state, plus a considerable amount of ran- 
dom noise, and ends at a state that shows features of 
both dipoles and bars. A tentative interpretation of this 
evolution is that the nonlinearity is simply exhausted by 
the viscous decay before either evolution can be com- 
pleted. Once the amplitudes fall below amounts or in 
Fourier configurations at which the activity is effectively 
nonlinear, the pattern in place is "frozen" in its topology, 
and can only slowly decay. The scatter plot has features 
of both the point and patch predictions. 

b. Several other runs demonstrated odd features. 
Metastable states were found that would persist for a 
long time as a consequence of evolution that, though dis- 
ordered, might be thought to be less than totally tur- 
bulent. An example appears in Figs. |l]. In this run 
(\/v = 12, 500), we began with a 16-pole solution of the 
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FIG. 21: A second run, whose late-time state is outside the patch/point classifications, that seems to have reached some 
metastable late-time state. Here, a large area of zero vorticity was created in the initial conditions by removing, asymmetrically, 
four pieces of a 16-pole patch solution. 
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FIG. 22: The ui — ip scatter plot of the late-time state achieved 
in Figs. |2l], suggesting two independent co-existing sinh- 
Poisson states which have developed asymmetrically. This 
is for the evolution shown in Figs. gjj . 



3-level patch equation, removed some of the squares of 
vorticity, and used the rest as an initial condition. Most 
of these went to the dipolar states, but one of the un- 
typical ones reached a vaguely quadrupolar one, in the 
last panel of Figs. ^l|. Fig. p2| is a w — ip scatter plot, 
showing a bilinear form, as if two locally "most probable" 
states had formed into a quasi-equilibrium that was slow 
to break up and decay toward anything globally "most 
probable." Figs. |23| show the evolution of the global 
quantities for this run. 

Several runs (1/z/ = 10,000) were carried out using 
this metastable quasi-quadrupole as the basis for initial 
conditions, plus a significant amount of added random 
noise. The evolution then was typically that the system, 
after having lingered awhile, evolved into the dipole con- 
figuration. 
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FIG. 23: Time evolution of the global quanties for the run 
shown in Figs. |l] and ||. Note that the mean wave number 
never reaches above 2.5. 



IV. CONCLUSIONS AND DISCUSSION 

We have set out to test the relevance of the predic- 
tions of Eqs. (||) and (||) and their respective vorticity 
discretizations to the numerically-determined, long-time 
states of a 2D NS fluid with Reynolds numbers of sev- 
eral thousand subject to doubly-periodic boundary con- 
ditions. Eq. (|^) has been derived by modeling the vor- 
ticity distribution as a mean-field limit of equal delta- 
function point vortices, while Eq. (^) models the vorticity 
distribution as made up of flat mutually exclusive patches 
of an area whose choice is arbitrary. Both equations have 
an infinite number of solutions, characterized by different 
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topologies, but for fixed vorticity fluxes and total energy, 
all but one are only local maxima and there seems to 
be always one uniquely denned maximum-entropy pre- 
diction. As seen in Figs. || through || it is sometimes a 
matter of painstaking analysis to determine which state 
this is, however. We did computational runs of two basic 
types: runs with broad-band, totally disordered initial 
conditions of the type previously investigated, and runs 
originating from vorticity distributions that consisted of 
large areas of nearly flat vorticity patches plus random 
noise to break possible unwanted symmetries. In the for- 
mer case, there were no surprises: the previously found 
dipolar late-time states inevitably resulted. This was 
not the case for the second kind of initial conditions, 
however. They did not generate broad-band turbulence 
in which energy was shared widely among many Fourier 
modes, but did exhibit considerable nonlinear activity in 
which the energy remained primarily in the lower parts 
of Fourier space. It is perhaps a semantic quibble as to 
whether the evolution should be called fully turbulent. 
This second class of initial conditions, in any case, ex- 
hibited a more diverse range of possible behavior than 
the broad-band initializations did, and often came close 
to a late-time state that could be identified as one of the 
solutions of Eq. (||). Most interesting of those was the 
one-dimensional "bar" state exemplified by Figs. || and 
|l3| ; in the former case, it is attained with the help of 
added initial random noise, and in the latter, without it. 
As Figs. H through^ illustrate, it may or may not be con- 
sidered to be the most-probable patch state, depending 
upon the choice of the patch size A/M, which seems to be 
arbitrary. In any case, it is not the most probable point 
state predicted from Eq. (Q). In the finite-size patch the- 
ory leading to Eq. , the arguments for it proceed from 
"coarse graining" the Euler equation behavior and pre- 
suppose a minimum observational scale below which fine 
spatial structure cannot be resolved. In a well-resolved 
NS computation (presumably including all of the ones 
reported here), there is no such scale, and observation of 
all scales participating is feasible. Thus that the patch 
formulation of Eq. (|6|), with a large enough value of A, 
can predict a radically different topology than it does 
with a smaller A, and then substantiate that prediction 
in a dissipative NS computation, is another of the accu- 
mulating puzzles which remain to be deciphered. In our 
opinion, there can be said to be an interesting 2D NS 
regime to be explored that has opened up in these inves- 
tigations in which nonlinear evolution, perhaps not fully 
turbulent in the conventional sense, nevertheless leads 
from one state that can be identified to another lami- 
nar, late-time state that is quantitatively more probable 
than its initial conditions: an essentially thermodynamic 
behavior. The conditions for assigning this probability 
remain incompletely defined, and seem to us worthy of 
further numerical investigation. 
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APPENDIX A: STATISTICAL MECHANICS IN 
2D TURBULENCE 

In this Appendix, we attempt to review the most el- 
ementary calculations of the "most probable" vorticity 
distribution, inside a periodic square, that is compatible 
with a fixed value of the total kinetic energy and fixed 
(equal) fluxes of positive and negative vorticity. We do 
this for the case of Lynden-Bell statistics, which can then 
be specialized to Boltzmann statistics as a special case. 
It should be emphasized that in both pictures, vorticity 
is regarded as a pointwise conserved quantity, so that 
viscosity is not involved. 

Consider the problem of rearranging a continuous pe- 
riodic vorticity distribution w(x, y; t — 0) to constitute a 
"more probable" arrangement of the vorticity elements. 
This can only be approached through some kind of dis- 
cretization. We consider coarse graining the vorticity dis- 
tribution in the following sense. 

Cut the periodic box into many cells of equal volume; 
i will label the cells. The cells are of area A. The coarse- 
grained vorticity distribution will be defined by the num- 
bers 

= x / / u(x,y)dxdy. (Al) 

A J J cell i 

We now represent the vorticity LUi inside the zth cell 
in terms of smaller sub-units. These can clearly either 
have finite areas ("patches") or zero areas ("points"). 
We should be able to take the "patch" formulation and 
achieve the "point" formulation by suitable limiting pro- 
cesses, so we start with the patches. The ultimate ques- 
tion will be which, if either, satisfactorily represents 
the long-time state of the Navier-stokes evolution of the 
u>(x, y;t = 0) at high Reynolds number. It is an aca- 
demic and not-very-interesting question as to which one 
better predicts the long-time Euler equation evolution of 
u>(x, y; t = 0), since we can never compute that anyway. 

Note that if uj(x, y;t — 0) is an analytic function and 
we look at the set of points x, y at which u> takes on 
any value u>q, say, that set is always a set of measure 
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zero i.e., zero area. So arriving at any finite number 

of "levels" will always implicity assume that some kind 
of "coarse-graining," as in Eq. (Al), has already been 
done. 

We now cut each cell up into equal areas (they do not 
have to have any particular shape), each one the size 
of a "patch." We say Mj is the maximum number of 
patches allowed in the ith cell. Without apparent loss of 
generality, we can take all the patches as of equal area, 
but we will carry along the index i on Mj. There is 
no obvious theoretical basis for choosing the size of a 
"patch," A/Mi, except that it is far smaller than A, so 
Mi > 1. 

We assume there are different "strengths" of the 
patches, though this will later come to seem perhaps un- 
necessary. Let the superscript j represent the "species," 
or strength of vortices of type j, which will be called 
K 3 . K 3 is the integral of the (uniform) vorticity of the 
patch over its area. Let j = 0, l,2,3...g, where K° = 
(an "empty" site, with no patch in it), and q different 
strengths represented, which can be positive or negative, 
but not zero. For convenience in notation, we will let 
K 3 = Kj, hereafter, to avoid confusion with exponents. 

Giving a set of occupation numbers Nf will give a dis- 
crctizcd microscopic representation of all the cof. 



Obviously, there is no unique way to do this. 
Note that 

N 3 < Mi and ^ n( = Mi. 

3 = 1 j=o 



(A2) 



Also Y^i — N 3 — const, for any rearrangement that 
does not create or destroy patches. will always run 
over all the cells.) 

For the "probability" of any total set of occupation 
numbers N 3 , we adopt the expression 



W = Const. \l JJ 



(A3) 



i i=i(JV?!)(Mi-£tfi)l 
i=i 



Notice that the second of the constraints ( |A2| ) has been 
built in and the N 3 may be varied independently for j — 
1,2, 3,4... q, so long as we demand that in the variation 
N 3 — > N 3 + 5N 3 

N? = Const. = N j ,j = 1, 2, 3, ...q, 



that 



£<Wf =0. 



(A4) 



The last term in the denominator of ( |A3| ) comes from 
the empty sites or "holes" in which no finite-vorticity 
patch resides. 



Assuming that all the N 3 ^> 1 and Mi ^> 1, we have 
for the entropy, using Stirling's approximation, 



S = lnW~Co72st.-J2j2(N 3 \nN 3 ~N 3 ) 

i j = l 

i ;=i i=i 

i 

-{Mi - JVj)} + (negligible terms). 



The discretely represented energy is 

E =\miY, K i N i^ K ^ 



(A5) 



(A6) 



i k j.l 



The j = and I = terms will give nothing, since Kg = 
0. Note that i/in- = ip^i is just the interaction energy of 
two unit-strength vortices in cells i and k, and that f/'ife 
is species-independent. 
We may also write 



where 



(A7) 



(A8) 



l.k 



ipi is the discretized stream function. 
When N 3 -> N 3 + SN 3 



8E = Y,Y, K i^ m i- 



(A9) 



For 5S, we have 

9 



ss = E E SN ? ln N i + E E ln ( M > - E N l) SN i ■ 

i j=l i j = l 1=1 

(A10) 



To maximize S under the constraints (A4) and SE = 0, 
using lagrange multipliers, we set 

i 

5S + f3SE + J2 E &3 ' SN i = °> 

i j=l 

where the Lagrange multipliers are (3 and &j , and the SN 3 
are now regarded as independent, so their coefficients 
may be set equal to for all i and j = 1,2,3, ...q. We get 



- ln Nj + &i + pKiipi + ln(M ( - ^ AT?) = 0, (All) 

we let f3 = —(3, ctj — otj to conform with earlier notation 
(they are still undetermined, anyway). 
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■j 1.2.:'....,/. (A12) 



Exponentiating 

Mi -Y,N\ 
1=1 

i 

But Mi — Nl — = number of empty sites in cell i, 
i=i 

always. 
So 



where lie inside the ith cell. In this formulation, 

the constraints are 



= E, 

i 

t j 



(A22) 
(A23) 



N j = N°e aj ~ ,3Kj ^ i 



and 



and so 



3=0 



3=0 



N? = Mi 



1=0 



j = 1,2,... q. 



(A13) 



(A14) 



(A15) 



We need a tractable expression for S, in order to com- 
pare different entropies for different sets of Nf for the 
same E and NK We may rewrite Eq. (|A3|) as 



i j=0 JV i ■ 



(A24) 



so that 



The same expression will work with j = if we take 
cko = 0. So 



N? = Mi 



1 



(A16) 



J2 eai-PKrft 
1=0 



and Eq. (A15) c an be e xten ded to the case j = 0, also. 

The equation ( A15| )-( A16 ) give the most probable set 
of occupation numbers 2Vj? for fixed E and iV J ' . If we take 
the Nf from Eqs. flAl5|) and ( A16 ), we have 



W=Y,Nt,j = 1,2,3. 



(A17) 



and 



5 = Const. - (-ZV? In N? -N?) + (negligible terms) , 

i 3 =0 

(A25) 

or 

S = Const. ~J2J2 N3 i ( a J ~ $ K 3^ + ln N i) 

i j = l 

-^ATPlrxATO + ^iVP, (A26) 

£ i 

5 = Const. - ^ ajiV J + 2/3S 

3 

- £(ln AT°)(£ ivf + JVP) + £ AT°. (A27) 

i j'=l i 

Notes that £| =1 iV/ + TV? = Mj, and =total un- 

occupied sites, and both of these are just fixed constants. 
So up to unimportant additive constants, the entropy of 
the stationary states is 



E =\Y.Y, K o N i^ K i N l (A18) S eq . = -^2a j N^ + 2l3E-J2M i lnN^ (A28) 



as q + 1 equations which determine the 9 + 1 unknowns 
atj and /3, j = 1, 2, 3 . . . g, as functions of E and the iV J . 

If we want associate a "partial vorticity" with each 
species, we may write 



Mi , 
1 = ~K * 



where 



Zi = £e°"-^S (oo=0). 



(A19) 



(A20) 



1=0 



Now the physical vorticity associated with the most 
probable state is 



The first two terms of Eq. ( |A28| ) are familiar from 
the "point" theory. The last one is new, and is a conse- 
quence of the finite area of the patches. If the "patch" 
area shrinks to zero, the last term becomes an infinite 
additive constant, apparently, and the "point" entropy is 
recovered. 

If we assume all the Mi are equal (there seems to be no 
obvious reason for not doing so), the last term simplifies: 

- Y M i ln N i =-J2(M\n M) + M In Zi ■ 

i i i 

The first term is constant, and the second term has the 
integral representation 



U)\ dxdy, 



(A21) 



M Y ln z i 



M 
~A 



whole box 



celli 



3=0 

(A29) 
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Note that M/A is still entirely up to us to choose. 

There is still far more freedom than we need to repre- 
sent any initial LOi. We can still do it, for example, by 
using only three value ol K: namely, Kj = —1, 0, or, +1. 
Any coarse-grained "level" can be achieved, that is, with 
only one kind of positive patch, one kind of negative one, 
and a zero or "empty" site. Then we have 

S = -a+N + - a_iV_ + 2f3E + M^lnZ,, (A30) 

i 

or in the continuous representation, 
M 



S = 



dxdy ln[l + e a+ -^ + e a - +M '} 



(box ) 

-a + N + - a_7V_ + 2/3E. 



(A31) 



If we now demand that ip obey Possion's equation, 
passing to the limit of arbitrarily many patches of ar- 
bitrarily small area, 



V 2 ip = -ui= — - 
A 



e a+-/3i/> _ e a_+/3i/> 



\ _|_ gQ + -/3l/> _|_ gQ_+/3V 



(A32) 



This would be the analogue of the "sinh-Poisson" or 
"tanh-Poisson" equation discussed later. 

It is simple, computationally, to consider initial states 
for which there is, for every u)i, another level with the 
value — Wj. This will in fact guarantee J J/ box \ Ludxdy = 0. 
It also makes it plausible to assume +, — symmetry, so 
that a + = at- . Then we would have 



M 
A" 



2sinh(/3V>) 
a + 2 cosh(/?V) 



(A33) 



Suppose we want to let the patch area — > 0, now. This 
would mean M/A — > oo . The only way the right hand 
side could stay finite and have a finite flux, J J \uj\dxdy < 
oo, would be for e~ Q — > oo, and do it in such a way that 
(M/A) / e~ a stayed finite: 



M 



A 2 

— = Const. 
Ae~ a 2 



So 



V 2 ^ = A 2 sinh(/3V), 



(A34) 



and the sinh-Poisson equation is recovered. 

It appears that for the symmetric periodic case, we can 
also get by with no unoccupied sites, that is N® = 0, Vi. 
Every Ui could be made up of just the right number of 
K = +1 and K = — 1 patches to give all the desired 
Ui, including zero. To treat this case, we can re-start the 
whole procedure, with only one species independent (K = 
+1, say), so that N~ = Mi — Nf, always. 

Then 

S = lnW^-J2(Nt^Nt-Nt) 

i 

- J2[(Mi - N+) ln(Mi - N+) - (Mi - N+)\, 



and the analogue of Eq. ( A12| ) is 

i {M l -N+), 



N+ = e Q -^* 



iV 4 



with 



AT," 



l+e°-* ~ e |(a-») + e -5(«-Wi)' 



This gives an overall vorticity of 

Mi smh(i(a-/3Vi)) 



A cosh(i (a -/%))' 



(A35) 

In this case, E is invariant to letting ipi — > tpi + const., 
so we can always choose the constant so that —(3 x 
Const. + a = 0, leaving 

Ui = — — tanh-^i. 

Since (3 is as yet undetermined, we might as well let (3 — > 
2/3, M,/A = A 2 , and we have 



V 2 ^ = -w = A 2 tanh(/?V'). 



(A36) 



Eq. ( A34) corresponds to the limi ting case in which the 
"patches" have no area, and Eq. (A36) to the limiting 
case in which they have all the area, and there is no 
empty space. Both apply to the symmetric case in which 
for every intial u>i, there is somewhere another with — u>i 
as its value. 

The most general, completely unrestricted, case has, 
as its analogue of sinh-Poisson, 



2^ a" J' "« 



(A37) 



1=0 



Comparing Eqs. ( |A37| ), ( |A3^ ), ( ggg ), and ( |A3^ ), it 
seems clear that every choice of the M/A and the pos- 
sible Kj will lend to a different equation and a different 
if). There seem to be no compelling arguments as to why 
one might be superior to the other in predicting long- 
time Navier-Stokes behavoir, and no point in speculating 
which will do a better job of predicting long-term Eulcr 
equation evolution, since we can never know that. It is 
perhaps worth remarking that if we want to shrink A 
so much that there i s no approximation involved in let- 
ting the Lot of Eq. (Al) represent an analytic function 
w(x,y,0), then the "points" are the only pos sibi lity. Ev- 
ery "level" approximation of the type Eq. ( |Al[ ), with i 
bounded from above is already a function significantly 
different from u(x, y, 0). For example, it has an infinite 
palinstrophy, and thus an infinite intial enstrophy decay 
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rate, because of the sharp discontinuities between the 
cells. 

A possible stategy is to search for intial conditions 
where the presence of finite-area conserved "patches" 
woul d see m to have the most important consequ ence s: 
Eq. ( A33 ) with positive a. For example, in Eq. ( |A33| ), 
which presumably can be solved numerically, choose an 
"a" and a "/? " wh ich will give a relatively flat scatter 
plot (like Fig. |l(b)l ). 

Pick an M / A that is large enough that tp seems slowly 
varying on the scale ~ y/~A, but M ^> 1. This will give 
a function uj(x,y) which is presumably well-represented 
by patches, if anything can be. Do not pick the maxi- 
mum entropy state but something like a quadrupole or 
higher multipole which can be expected to break up and 
generate turbulence when used as an initial condition for 
Navier-Stokes code at small v. 

There may have to be a small random perturbation on 
those initial conditions to get it to go turbulent as soon as 
possible. Notice that choosing a, j3 and M/A in advance 
provides computable expressions for energy and vorticity 
fluxes. The first question to ask is, what "final" state 
does the system relax to for large times - what does its 
scatter plot look like? 

The second question is, then, using those expressions 
for energy and vorticity flux in the sinh-Poisson code, 
what does the final maximum-entropy state look like - 
its scatter-plot in particular? Is there any sense in which 
the initially flat areas of vorticity flux are still there? 
The same program coul d be carried out starting from 
the more restrictive Eq. ( A36| ). 

A really convincing, but perhaps over-ambitious plan 
would be to tak e the a, (3 and vorticity flux from the 
solution to Eq. ( |A33| ), and use them as input for gener- 
ating a sinh-Poisson solution, comparing its scatter plot 
with the one resulting dynamically. 



APPENDIX B: SOLVING THE POISSON 
EQUATION 

The sinh-Poisson equation (Eq. ([)[)) has been solved by 
a lot of people numerically and analytically (see Ref. |2(J 
and references therein). Probably the most direct way to 
do it is to put the equation into spectral space and do 
the iteration there: 



\2 - 

(*»+i)k = 7r^(sinh*„) k . 
k 



(Bl) 



We show three solutions in Fig. g. In addition to the 
solutions exhibited, there are an infinite number more, 
characterized by more and more maxima and minima. 
They are obtained by startin g w ith a trial function on 
the right hand side of Eq. (Bl) that has the desired 
number of final maxima and minima, as noted by several 
authors: McDonald @, Book et al. @, Lund gren and 
Pointin S. Typically, the solutions with many maxima 



and minima are lower entropy solutions, for fixed energy 
and vorticity fluxes, and so we do not discuss them in 
any detail. 

The spectral method loses it advantage when we try to 
use it to solving Eq. ([l5]), because we might have some 
step function solutions. Here we adopt the numerical 
method used by McDonald ||. 

The first step is to get the nontrivial solution for Eq. 
@J|) with * = on a square boundary of [0,tt] x [0, w]. 
Starting from a trial solution W(x, y), we get the solution 
of v(x, y) from 



V 2 v + vf'(W) =R(W,W), 
subject to v — on the boundary. Here 



(A,B) = / / 
J o Jo 



/(*) = 



A(x,y)B(x,y)dxdy, 



A 2 sinh * 



g + cosh ' 



R(x,y) 



V 2 W + f(W) 



(W,W) 

Then we can correct the trial solution: 

v(x,y)(W,W) 



w ->w 



2{v, W) - (W, W) ' 
until we get a sufficient accurate solution. 



(B2) 



(B3) 



(B4) 



(B5) 



(B6) 
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FIG. 24: A periodic quadmpolar solution may be used to gen- 
erate a dipolar one, depending upon the spatial sub-volume 
chosen, as shown. The value of A 2 in Eq. ( |l5| ) will changed 
correspondingly. The procedure is somewhat involved. 

Perhaps the most important change we made to Mc- 
Donald's methods H is using the quadruple precision 
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instead of double precision in the calculation (and this 
is not a luxury on mordern computers). The accuracy 
problem mentioned by McDonald in this method is over- 
come. As a result, if we define the root-mean-square error 
(RMS): 



where 



and 



e = (£/7£<? 2 ) 



1/2 



h3 i,3 



f = A^(i,j) + X 2 S mh^(i,j)) 



g= \AV(i,j)\+ |A 2 sinh(*(z,j))| 



Our solutions can at least reach an accuracy of 10 7 . 



We may construct higher-order multipole solutions by 
putting monopolc solutions side by side with alternat- 
ing signs, and achieve doubly-periodic boundary condi- 
tions in this way, in a "checkerboard" solution, so that 
dipole (Figs. g4|) and higher-multipole solutions can be 
constructed. There is some loss of accuracy in this pro- 
cedure, and we have made sure that the accuracy of all 
solutions is of the order of 10 -6 (RMS) in our calcula- 
tions. 
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